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ABSTRACT 

We  can  create  a  richer  and  more  neurophysiologically  realistic  model  of 
neural  activity  in  the  brain  by  developing  a  model  of  neural-dendritic  coupling, 
one  which  expressly  accounts  for  the  way  in  which  the  many  afferent  connec¬ 
tions  into  the  neural  body  influence  the  somatic  membrane  potential.  Such  a 
model  would  begin  to  fill  the  need  within  the  Artificial  Neural  Network  commun¬ 
ity  for  neural  models  which  go  beyond  the  current  "weighted  sum"  paradigm  for 
artificial  neuron  connectivity.  Although  such  models  have  use  in  engineering 
applications,  there  are  many  aspects  of  biological  neural-dendritic  organization 
which  could  enrich  artificial  neural  networks.  Moving  from  simple  "axonal" 
connection  weight  neural  models  to  neural -dendritic  models  with  a  richer  struc¬ 
ture  will  allow  investigation  of  both  events  at  the  neural  level  (e.g.  inter-spike 
interval  histograms  and  stochastic  resonance)  and  also  potentially  at  the  neural 
systems  level.  This  will  also  introduce  the  possibility  of  introducing  cross-scale 
interactions  into  artificial  neural  systems. 


I.  INTRODUCTION 

The  effect  of  neural -dendritic  interactions  has  so  far  been  only  weakly  probed  in  the  realm 
of  artificial  neural  networks  and  neural  modeling.  Traditional  Artificial  Neural  Network  (ANN) 
models  of  interacting  neurons  use  a  simple  description  of  neural  connectivity.  In  such  simple 
models,  communication  between  neurons  is  afforded  by  artificial  axons,  and  the  "strength"  of  a 
given  neuron-to-neuron  connection  is  given  as  a  "connection  weight”  between  the  two  neurons. 
Synaptic  plasticity  is  viewed  in  terms  of  the  modifiability  of  connection  weight  ("axonal") 
strengths  (for  a  review  of  ANNs,  see  Maren,  Harston,  and  Pap,  1990). 

A  more  neurophysiologically  rcalisitic  -  and  interesting  -  basis  for  modeling  neural  systems 
would  take  into  account  the  nature  of  neural -dendritic  connectivity.  Such  considerations  for 
neural  systems  modeling  were  advanced  as  early  as  1958  by  John  von  Neumann.  In  'The  Com¬ 
puter  and  the  Brain.'  von  Neumann  wrote: 

..However,  the  more  frequent  situation  is  that  the  body  of  a  neuron  has  synapses  with  axons 
of  many  other  neurons.  It  even  appears  that,  occasionally,  several  axons  from  one  neuron 
form  synapses  on  another.  Thus  the  possible  stimulators  are  many,  and  the  patterns  ot 
stimulation  that  may  be  effective  have  more  complicated  definitions  than  the  simple  "and 
and  "or"  schemes  described  ...  It  may  well  he  that  certain  nerve  pulse  combinations  will 


second-order  chemically  gated  synapses  have  longer  time  constants  and  cannot  he  treated  In  imr 
model. 

in  order  to  characterize  processes  at  the  synaptic  level,  tc  .  to  characterize  a  dciidntn 
volume  activation  in  terms  of  its  relationships  to  inputs  and  to  arm  anon  deca>  vo  -e>  .  ■  l  •  n 

\k  the  iotal  post-synaptic  activation  of  a  small  dendritic  volume,  which  can  include  several  den 
dritic  spines  and  their  afferent  connections.  The  dendritic  volume  is  rough!}  equivalent  to  the 
"dendron"  introduced  by  Eccles  (19641.  Pre-svnapuc  chemically-gated  input  to  the  dendntu 
volume,  resulting  from  action  potentials  at  other  neurons,  contnbutes  to  the  total  dcndntK 
volume  activation.  We  think  of  this  input  as  being  predominantly  due  to  chemically-gated  synap 
tic  transmissions  which  have  time -constants  in  the  l  -5ms.  range  (Kandcl  et.  a!  lh‘H  pa.  124  ;  ']  o 
a  great  extent,  the  elicited  post-synaplic-potennal  CPSPi  response  to  aflcretu  spiking  is  linear 
However,  the  response  is  bounded  from  above  in  die  case  of  excitatory  inputs  and  from  below 
when  the  inputs  arc  inhibitory.  This  is  due  to  probabilistic  release  of  quantized  packets  of  neuro 
transmitter  into  the  synaptic  cleft  (Eccies,  1981}  and  to  a  finite  number  of  channels  in  the  post- 
synaptic  terminus.  Further,  experimental  measurements  of  excitatory  and  inhibitors  PSPs  show 
that  both  are  maximally  within  several  tens  of  millivolts  of  the  membrane  resting  potential  [Kan- 
del  et.  al.  1991).  This  also  indicates  bounding  of  the  PSPs.  Bounding  also  comes  about  through 
scalingof  the  PSP  which  occurs  in  the  short  distance  the  PSP  travels  from  the  synaptic  site,  at  the 
spine  head,  to  the  spine  base;  simulations  indicate  that  the  potential  may  drop  by  as  much  as  959? 
as  a  result  [Segev  et.  al.  1992).  Consideration  of  these  combined  factors  will  allow  us  to  model 
the  temporal  derivative  of  the  activation  or  total  polarization  of  a  small  dendritic  volume  as  an 
activation  decay  term  plus  input  terms.  The  inputs  will  include  the  weighted  and  bounded  effects 
of  presynaptic  signals  from  axonal  and/or  dendritic  connections,  a  weak  (low  amplitude  and  fre¬ 
quency)  periodic  driving  force,  and  noise.  As  we  shall  see  later,  the  form  for  the  temporal  deriva¬ 
tive  of  the  activation  in  a  synaptic  volume  will  be  the  same  as  for  the  neural  body  itself 

The  activation  potential  which  originates  in  a  local  dendritic  volume  and  moves  through  the 
dendritic  arborization  towards  the  soma  undergoes  two  substantial  changes  [Segev  et.  al.  1992], 
Both  changes  influence  the  model  that  we  construct  for  time-dependent  activation  at  the  soma. 
First,  there  is  a  marked  attenuation  of  the  maximal  signal  amplitude,  and  second,  there  is  a  draw  ¬ 
ing  out  of  the  signal  waveform.  Simulations  by  Rail  and  Rinzel  [1973],  Rinzel  and  Rail  [1974] 
and  Segec  et.  al.  [1992]  indicate  a  strong  attenuation  if  the  synaptic  PSP  as  it  passes  from  the  ori¬ 
ginating  dendritic  volume  towards  the  soma.  The  final  maximal  signal  amplitude  reaching  the 
soma  can  be  less  than  1%  of  the  original  maximal  amplitude,  even  though  the  stretched  signal 
duration  indicates  that  a  substantial  portion  of  the  signal  is  preserved  via  temporal  integration. 
This  leads  us  to  model  the  neuronal  (i.e.  soma)  input  not  as  the  direct  value  of  the  synaptic  activa¬ 
tion  but  rather  as  a  bounded  function  of  it.  This  makes  plausible  the  use  of  a  transfer  function 
such  as  a  sigmoid  or  hyperbolic  tangent,  which  is  so  ubiquitous  in  neural  modeling  It  is  impor¬ 
tant  to  point  out  that  the  decoupling  (to  be  described  in  the  following  section)  of  the  neuro- 
dendritic  stochastic  differential  equations  via  the  adiabatic  theory  will  be  tantamount  to  a  quasi- 
linearization”  of  the  dendritic  dynamics.  Another  observation  arising  out  of  these  simulations  is 
that  the  PSP  signal  waveform  is  greatly  stretched  as  it  travels  towards  the  soma  Specifically,  an 
input  that  may  occur  in  less  than  10ma.  at  the  synaptic  site  may  have  an  influence  persisting  over 
lOOms.  at  the  soma.  This  is  particularly  true  if  the  synaptic  site  is  far  from  the  soma  We  can 
summarize  this  by  noting  that  "...the  dendritic  tree  behaves  as  a  substantial  dcla\  line  for  the 
synaptic  inputs  [and  that]  distal  inputs  are  subject  to  significantly  longer  neural  delays  than  proxi¬ 
mal  inputs..."  [Segev  et.  al.  1992],  The  notion  of  different  distributions  of  synaptic  input  has  been 
stated  and  experimentally  observed  as  early  as  the  1960s  |e.g.  Rail  1970J.  In  am  given  neuron, 
due  to  the  branching  nature  of  the  dendritic  tree,  there  will  be  far  more  distal  input--  than  proximal 
inputs.  Thus,  to  a  first  approximation,  we  can  model  the  synaptic  inputs  under  the  assumption  that 
the  time -constants  of  events  are  "stretched”  as  they  move  through  the  dendntu  passage  to  the 
soma  This  allows  us  to  make  an  adiabatic  approximation  in  treating  the  actw.i;:  •:  equations  or 
the  dendritic  volumes  and  the  neuron. 

We  turn  our  attention  now  to  the  membrane  potential  at  the  neurai  sen:. a  '  .-.ve.  hillock, 
and  note  that  again  we  ma\  think  of  the  temporal  derivative  ('I  the  a.  mao--  ”  'em 


hack  to  a  (zero-level)  resting  state  A  more  effective  model  of  neural  activation  recjutres  th.i-  the 
inputs  to  the  neural  system  he  modeled  tn  a  continuous  manner,  consistent  vwtli  our  understand 
me  of  the  ".stretched"  arrival  n!  dendritic  activations  at  the  neural  soma,  as  has  been  described  h\ 
Rail  and  Rinzcl  1 19731-  R in/el  Ac  Rail  1 1974 j.  Holmes  A  Rail  i  l‘>92|.  and  Seges  et  al  ;  !<»>:;  hi 
a  later  work,  Stein  et  al  1 1474;  identified  the  activation  oi  the  neuron  as  the  rate  at  which  nerve 
impulses  would  be  ft  red  While  this  would  allow  an  intuitive  connection  between  the  neural 
model  description  and  "states",  which  could  be  roughly  described  in  terms  of  the  average  fre¬ 
quency  of  firing  of  an  action  potential,  it  dees  not  serve  when  the  model  needs  to  be  connected 
with  a  more  explicit  description  of  neural  dynamics,  such  as  the  generation  of  intersptke  interval 
histograms.  We  will  return  to  this  point  later  in  tins  work.  Subsequent  investigators  have 
developed  diffusion  models  in  which  the  discontinuities  in  u,(n  have  been  smoothed  out  (see  e  g 
Tuck  we 11.  1980  and  references  cited  therein).  While  more  mathematically  tractable,  such  models 
again  treat  the  arrival  of  excitatory  and  inhibitory  signals  to  the  neural  soma  as  discrete  rather 
than  continuous  processes. 

For  the  current  work,  a  more  precise  identification  of  the  neural  membrane  potential  is 
needed;  one  which  distinguishes  it  from  the  potentials  within  the  dendritic  network  afferent  to  the 
soma,  and  one  which  decouples  die  continuous-time  activation  changes  in  the  soma  that  occur  in 
response  to  input  and  activation  decay  from  the  action  potential,  whose  abrupt  nature  may  be 
viewed  as  a  reset  mechanism  of  an  otherwise  continuous  process.  In  our  work,  the  variable  uAt ) 
refers  specifically  to  the  membrane  potential  at  the  trigger  zone  in  the  neural  soma.  This  is 
because  the  interesting  dynamics  of  the  soma  are  generated  at  the  trigger  zone,  which  has  a  lower 
threshold  than  the  rest  of  the  soma.  However,  changes  in  membrane  potential  at  the  trigger  zone 
propagate  rapidly  both  throughout  the  soma  itself,  and  down  the  axon  as  an  action  potential.  We 
regard  the  brief  depolarization  and  ensusing  hyperpolarization  of  the  action  potential  (and  its 
corollary  within  the  neural  soma  itself)  as  a  reset  mechanism  whose  details  are  not  addressed  in 
this  work.  Our  model  does,  however,  address  the  continuous  changes  in  membrane  potential  due 
to  dendritic  connectivity,  activation  decay,  and  other  factors. 

The  disparity  in  the  neural  and  dendritic  time-constants  discussed  above  is  incorporated  into 
our  model  via  the  constraint, 

/?,«/?,  (i >1) .  (2) 

This  constraint  will  allow  us  to  confidently  utilize  the  slaving  principle  of  Haken  [  1977]  to  reduce 
the  system  (1).  This  reduction  is  earned  out  in  the  next  section.  Wt  also  consider,  via  numerical 
simulations,  the  range  of  validity  of  the  adiabatic  elimination  technique  as  well  as  the  bifurcation 
properties  and  collective  effects  that  appear  in  the  decoupled  neuron  dynamics  because  of  the 
interaction  with  the  dendritic  bath.  Finally,  we  discuss  switching  events  (between  the 
firing/quiescent  states  of  the  neuron)  that  arise  as  a  consequence  of  the  interplay  between  the 
noise  and  the  periodic  modulation...  stochastic  resonance. 

II.  SINGLE  NEURON  DYNAMICS 

The  separation  of  time-scales  embodied  in  the  inequality  (2)  permits  us  to  apply  the  slaving 
principle  [Haken  1977]  to  the  coupled  system  (1).  This  leads  to  a  closed  equation  for  the  soma 
activation  function  u{(i)  i.e.,  the  system  (1)  is  decoupled.  The  proceedure  is  outlined  in  the  fol¬ 
lowing  subsection  (refer  to  Schieve,  Bulsara  and  Davis  1991,  Buisara,  Maren  and  Schmera,  1992 
for  details)  in  which  we  also  introouce  the  effective  "potential  function"  corresponding  to  the  sin¬ 
gle  neuron  dynamics.  This  is  followed  by  an  analysis  of  the  bifurcation  properties  of  the  reduced 
model  (in  the  absence  of  the  deterministic  modulation  q  sinew).  The  remainder  of  this  paper 
examines  the  cooperative  effects  (together  with  their  potential  implications  in  neuroscience)  that 
arise  when  this  modulation  term  is  switched  on. 

Reduced  Neuron  Equation 

Equation  ( 1 1  represents  a  system  of  globally  coupled  nonlinear  stochastic  differential  equa 
lions,  subject  to  the  constrain;  (2i  on  the  equivalent  circuit  resistors  Rt.  The  N -dimensional 


where  o~  &  or  *  and  fa  i  is  Gaussian  delta-correlated  noise  having  zero  mean  an. 
and  we  have  introduced  the  (deterministic i  potential  (unction  tliat  allmds  u\  a 
mathematical  treamient  of  the  dynamics: 

l  u<  i  ^  an  '  p  In  coshu 

The  coefficients  m  the  above  dvnamtes  are  eiven  bv. 
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and  we  have  set  Gk  =  1-AU  Equation  (9)  (the  so-called  "reduced"  or  'effective'  neuron 
dynamics)  constitutes  the  starting  point  for  our  subsequent  analyses,  u ,  representing  the  voltage 
at  the  soma.  We  note  that  it  contains  (via  (11))  the  effects  of  the  coupling  between  the  cell  body 
and  each  of  the  dendritic  volumes.  Back-coupling  (characterized  by  the  coefficients  jt  ,  i  effects 
are  included  as  well  as  self-feedback  terms  (characterized  by  the  coefficients  J„).  Cross-coupling 
terms  (characterized  by  products  of  the  form  RlRkJlkJkl ;  /ta  >l)  between  the  dendritic  spaces  are 
also  present  but  are  neglected,  ft  is  important  to  point  out  that,  in  general,  the  expressions  ( 1  Ib.ci 
contain  higher  order  terms;  these  are  absent  because  of  the  truncation  (at  second  order)  of  the 
Taylor  expansion  of  the  functions  tanhu,*,,  consistent  with  our  assumption  of  small  deviations 
from  equilibrium  of  the  dendritic  "bath".  The  detailed  analysis  [Schieve,  Bulsara  and  Davis  1991 1 
leading  up  to  the  expressions  (11)  shows  these  terms  to  be  of  higher  order  in  o}R,i2C,.  We 
assume  that  the  dendritic  noise  terms  are  sufficiently  weak  so  that  we  always  have 


2C, 


<1  i  >1 


(12) 


so  that  the  deviations  u<  -u,  are  small.  Note  also  that  we  have  not  assumed  that  the  matrix  J  is 
symmetric,  in  contrast  with  existing  (Hopfield-type)  models.  In  fact,  no  further  assumptions 
beyond  (2)  and  (12)  (both  of  which  are  based  on  very  reasonable  physical  and  neurophysiological 
arguments)  as  well  as  the  adiabatic  assumption  of  very  low  modulation  frequency  co  need  be 
made  to  obtain  the  reduced  neuron  dynamics  (9).  In  the  next  subsection,  we  write  down  the 
"potential"  function  corresponding  to  the  reduced  neuron  dynamics  and  examine  its  stability  and 
bifurcation  properties.  We  also  present  numerical  simulation  results  that  elucidate  the  range  of 
validity  of  the  reduced  description  (9)  when  compared  to  the  full  dynamics  (1). 


Steady-State  Potential  Function;  Stability  and  Bifurcation  Properties 

It  is  easy  to  show,  via  linear  stability  analysis,  that  the  dynamics  in  (9)  are  globally  stable  as 
long  as  a>0.  For  this  case,  the  potential  function  U (u ,)  is  a  Liapounov  function  for  the  reduced 
dynamics  (9).  We  note  that  stability  does  not  depend  on  the  properties  of  the  coupling  matnx  J 
Assuming  a>0,  we  can  easily  show  that  the  potential  (10)  will  be  parabolic  (with  an  elliptic  fixed 
point  at  uj=0)  if  p<a  (this  includes  negative  as  well  as  positive  values  of  p).  For  p  >  a,  the  poten¬ 
tial  is  bimodal  (with  its  minima  located  at  c  =±(p/a)tanh(p/a),  and  a  hyperbolic  fixed  point  at  ,  mm 
and  cooperative  stochastic  effects  come  into  play.  The  transition  to  bimodaliu  (at  p-cn  is 
accompanied  by  a  pitchfork  bifurcation  in  the  most  probable  value  of  the  activation  u  with  the 
two  states  (attractors)  corresponding  roughly  to  the  quiescent  and  firing  states  ol  the  neuron.  The 
flow  (given  by  the  first  term  on  the  rhs  of  (9))  exhibits  the  characteristic  N-shaped  relationship 
known  to  exist  in  excitable  cells  (see  e.g.  Rinzel  and  Ermentrout  1989;  Abbott  and  Kepler  I  non. 
We  now  consider  the  effects  of  the  cell  body  coupling  to  multiple  dendritic  volumes  on  -  t:.r 
sition  Throughout  the  remainder  of  this  work  we  shall  assume,  for  simplicity,  that  the  rim-e  \  .tr¬ 
ances  in  the  elemental  dendritic  volumes  are  the  same  and.  further,  that  all  the  domino, 
constants  are  equal  Specifically,  we  set  or;  -•  ct?  ,R,  -R:  for  all  t  >  l  and  C  !  fora’'  w  • 


(ln)d 


We  now  digress  brief!)  to  consider  she  case  in  which  the  potential  is  monomodai  m  the 
absence  of  any  coupling  to  the  dendritic  bath  (this  can  be  achieved  h\  setting  /  K  <  !  i.  m  con- 
trast  to  the  case  discussed  in  the  preceding  paragraph  Then,  one  can  easils  calculate  the  value  of 
R :  tier  given  noise  variance  05  and  configuration  of  the  matrix  Ji  above  which  the  effector 
potcruial  is  bimodal.  Increasing  R,  leads  to  a  transition  to  bimodality  (occurring  at  [>  o  l »  nnl \  I01 
the  case  in  which  the  sum  £ ./ ;  is  positive  (keeping  in  mind  the  constraint  imposed  b\  the  ine¬ 
quality  (12)).  This  may  be  realized  by  imposing  the  same  sign  on  the  vast  majority  of  the  off- 
diagonal  elements./,,  and  7lS.  Increasing  the  noise  variance  a;  degrades  the  effect,  this  is  evident 
from  (lib).  U  is  apparent  that  the  coupling  to  the  dendritic  bath  may  actually  introduce  a  phase- 
transition-like  behavior  into  the  neuron  dynamics  Effects  such  as  this  couphn^-uuiu ced  bimodal¬ 
ity  are  a  hallmark  of  multiplicative  noise  (see  c.g.  Horsthcmke  and  Lcfcvcr  19X4 1  and  have  been 
examined  in  simpler  neuron  models  (for  both  additive  and  multiplicative  noises i  by  Bulsara,  Boss 
and  Jacobs  [1989],  and  Bulsara  and  Schieve  [1991],  The  opposite  effect  can  also  occur:  depend¬ 
ing  on  tite  magnitude  and  sign  of  each  element ./,,,  a  potential  that  is  bistable  in  the  absence  of  the 
bath  coupling,  can  be  rendered  monostablc  by  the  dendritic  field.  This  is  evident  from  (11b). 
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j  Fig  2.  Neuron  (cell  body)  probability  density 
1  computed  via  direct  integration  of  (1)  and 
i  the  reduced  equation  (9)  (data  points). 
\(R\,af,g2,q  ,  to, /VJ  =  (10,5,2,0.1,0.1,10). 

!  /?2  =  0.35  (solid  curve  and  filled  data  points) 
j  and  1.0  (dashed  curve  and  asterisk  data  points), 
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Before  concluding  this  subsection,  we  present  the  results  of  numerical  simulations  aimed  at 
demonstrating  the  validity  of  the  approximations  made  in  this  paper.  Figure  2  shows  the  probabil¬ 
ity  density  function  corresponding  to  the  "slow"  variable  u,.  The  solid  curves  (corresponding  to 
simulations  of  the  fully  coupled  system  (1))  have  been  obtained  via  direct  integration  on  an  HP- 
Apollo  425T  workstation.  We  consider  Av =10  elements  since  the  simulations  become  prohibi¬ 
tively  time-consuming  for  larger  N  values.  The  system  has  been  integrated  through  12,0(X) 
periods  of  the  deterministic  modulation  (after  allowing  the  transients  to  die  out)  using  a  stepsize 
of  0.015.  Noise  has  been  included  in  the  dynamics  via  the  Heun  algorithm  (see  e.g.  Greiner  et  al. 
!  988).  The  data  points  show  the  results  of  integrating  the  corresponding  effective  one-body  equa- 
: ion  i‘J>  using  the  same  routine  In  this  figure  we  take  ~  10  and  consider  the  case  R2~ 0.35  and 
1  0.  all  other  parameters  remaining  fixed  The  agreement  between  the  exact  and  reduced  dynam¬ 
ic-  i  -  seen  to  be  excellent  for  R2  =  0.35  Increasing  R2  to  1.0  still  yields  reasonably  good  quahta- 
uv  aen-rriH’nt  although  it  is  apparent  that  we  arc  very  close  to  the  boundary  at  which  (2)  ceases 
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An  approximate  expression  for  the  power  spectral  density  lor  a  one-dimensional  nonlinear 
stochastic  system  of  the  form  (9)  has  been  derived  by  McNamara  and  VViesenleki  ;  i(»vr 

7’(D>  =  1  1  -2Z  (r0'r  1  (8Zr„e2)  -*■ -5/  Tt(>(1£r 2)‘S(Q-oj  i 

:■  N(Q)  -r  S(Q)5(£2  -  to> ,  (  l(ii 

where  we  define  Z  =  (4r(;  +  Q’)-!  and  ^  =  6 cg{  is  a  perturbation  theory  expansion  parameter  S 
and  N  arc,  respectively,  the  signal  and  noise  powers.  The  above  expression  provides  a  good 
approximation  to  the  power  spectral  density  for  £  <  1  and  has  been  derived  under  an  adiabatic  (i  c 
low-frequency)  approximation  to  «  Ua\c )  that  is  somewhat  less  stringent  than  the  to  «  r(,  approx¬ 
imation  introduced  earlier.  The  SNR  is  then  obtained  from 

r  r  s  i 

SNR  =  1 0  log  <1  -5^-  +  N  ( o) )  > !  .  (17) 
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1  Fig  3.  SNR  computed  from  (16)  and  (17)  for 
i  (R ,  ,R  2 ,  q  ±co ,  Aw  ,N )  =_£] 0X>.6,0. 1 A 1 ,0.00 1 ,100). 
"V i.  =  1  =-Jn  =  1  =Jl-  Ju  =0 ,yuJ=  1  ,(i  >1). 

""J  Bottom  curve:  Ju  =0  (isolated  case).  Remaining 
curves:  o22  =  0  (top),  1  (middle),  and  2  (lower), 
j  Inset:  the  case  for  no  modulation  in  the  cell-body 
1  eqn.  in  the  system  (1). 
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Fig  4.  Peak  SNR  (normalized  lo  its  value 
:  for7u  =0  case)  vs.  for  (R ■  ,a  w,Aw,.v  i 
l  (10,0.1,0.1,0.001,10)  and  cr?  =  0  (top  curve). 

-I  Kmiddje  curve),  2  (bottom  curx'c). 
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Solid  curves:  Ju  =  1  =7,.  j?t  -  i  =.JJ 
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absence  of  the  bath,  as  a  consequence. 


IV.  STATISTICAL  ANALYSIS  OF  FIRING  EVENTS 

It  is  important  to  point  out  that  stochastic  .csonance  (as  characterize'.!  by  the  bell-shaped 
SNR  vs.  noise  variance  curve  of  figure  3)  has  not  yet  been  directly  observed  in  a  living  system 
(although  we  will  ennunciatc  a  possible  caveat  to  this  statement  at  the  end  of  the  following  sub¬ 
section).  The  existence  of  noise-induccd  switching  in  the  nervous  system  would  seem,  however, 
to  be  an  eminently  reasonable  assumption,  based  on  our  simple  model  of  the  neuron  as  a  noisy 
bistable  switching  element.  Certainly,  noise  is  ubiquitous  in  the  nervous  system;  hence  one 
might  expect  that  when  sensory  neurons  arc  periodically  stimulated,  the  time  intervals  between 
successive  firing  events  contain  sensory  information.  These  "reset"  or  "refractory”  events 
correspond  to  the  repolarization  of  the  cell  membrane.  In  fact,  this  has  been  well-known  to  neu¬ 
rophysiologists  for  many  years.  In  neurophysiological  experiments  it  is  common  to  assemble  an 
ensemble  of  firing  events  and  fit  a  histogram  to  the  intervals  between  the  spikes.  An  example  of 
these  Inter-Spike-Interval  Histograms  (ISIHs),  which  are  quite  commonly  seen  in  the  neurophy¬ 
siological  literature,  is  shown  in  the  following  subsection.  In  this  subsection  we  also  demonstrate 
how  the  salient  features  of  the  experimentally  observed  ISIHs  can  be  easily  explained  by  our 
bistable  model;  in  fact,  we  shall  see  that  our  model  affords  the  simplest  possible  interpretation  of 
the  experimental  spike  train  data.  Throughout  the  rest  of  this  work  we  show  results  for  arbitrary 
values  of  the  constants  a, (3,6  in  (9),  i.e.,  we  do  not  consider  the  particular  details  of  the  quanti¬ 
ties  on  the  right-hand-sides  of  equations  (11). 

The  Inter-Spike-Interval  Histogram 

We  now  return  to  our  reduced  system  (9)  and  consider  only  the  time  intervals  of  the  transi¬ 
tions  between  the  potential  wells  (labeled  A  and  B)  while  ignoring  the  intrawell  motion  (recall 
that  the  potential  wells  correspond,  roughly,  to  the  firing  and  refractory  states  of  the  soma).  This 
is  tantamount  to  replacing  the  detailed  dynamics  contained  in  (9)  by  the  equivalent  "two-state 
dynamics"  depending  only  on  the  barrier  height  and  the  locations  of  the  elliptic  points  of  the 
potential  (10).  The  result  is  the  random  "telegraph  signal”*  (r)  (  £«,(()  in  the  notation  of  this 
paper)  depicted  in  figure  5. 


ABBA 


ABAB 


Fig  5.  Sample  output  from  two-state  device 
driven  by  noise  plus  weak  sinusoidal 
modulation.  Stable  states  at  ±*0  are  labelled 
A  and  B.  The  ABBA  and  ABAB  sequences 
shown,  arc  the  only  consecutive  time-interva 
sequences  available  to  the  two-state  system 
from  which  ISIHs  can  be  generated. 
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the  deterministic  modulation  ir.  :9>.  the  tv\ {»  sequences  of  figure  5  lead  to  different  ISIH-  'i  fie  top 
sequence,  referred  to  a>  the  ABBA  process,  measures  the  escape  time  from  well  B  It  k\uK  to  ,e, 
ISIH  with  modes  located  at  <><u!  integer  multiples  of  I  2.  T  heme  the  moduiatton  period  A  theor. 
describing  this  histogram  has  recently  been  developed  [Zhou.  Moss  and  June  PAM  ,  1  he  bottom 
sequence  (the  ABAB  process)  leads  to  ••  histogram  with  peaks  located  at  ail  micitc'  multiple'.  ■' 
7  This  is  the  sequence  commonly  oi"..rvcd  in  experiments  (and  the  one  that  we  concentrate  on 
through  the  remainder  of  this  livussion):  it  points  to  the  existence  ol  a  "reset  mechanism 
between  every  pair  of  spikes.  The  reset  events  arc  identified  w  ith  the  rcpolan/aiions  of  the  neuron 
membrane  that  occur  between  successive  upstrokes  of  the  action  potential  and  are  not  direct]', 
observable  in  ncuropnysiological  experiments,  in  figure  6  we  show  an  experimental  ISIH 
obtained  from  the  single  auditory  nerve  fiber  of  a  cal.  This  daia  should  be  compared  with  the 
.ABAB  ISIH  shown  in  figure  7.  The  sequence  in  this  figure  is  obtained  via  analog  simulation  ol 
(9)  with  the  potential  function  U  given  by  (10)  as  well  as  the  "standard  quarttc 

=  “t2  ^  j  “i*  Bins  potential  is  also  bistable).  The  sequence  of  peaks  in  the  ISIH  implies  a 

form  of  phase -locking  of  the  neuron  dynamics  to  the  stimulus.  Starling  from  its  quiescent  state, 
the  neuron  tries  to  fire  at  the  lirst  maximum  of  the  stimulus  cycle.  If  it  fails  to  do  so,  it  will  fire  at 
the  next  maximum  of  the  stimulus  (i.c.  after  a  complete  stimulus  cycle)  and  so  on.  with  a  finng 
event  corresponding  to  a  switch  between  the  two  states  of  the  potential  (10).  This  "statistical  skip¬ 
ping"  leads  to  the  sequence  of  peaks  in  the  ISIH.  Decreasing  the  noise  strength  (keeping  all  the 
other  parameters  fixed)  leads  to  more  peaks  in  the  histogram  since  skipping  becomes  more  likely 
Conversely,  increasing  the  noise  lends  to  concentrate  the  probability  into  the  first  few  peaks.  For 
vanishingly  small  stimulus  amplitude,  the  peaks  merge  into  a  Gamma  distribution  characterizing 
the  ISIH  for  the  spontaneous  case  [Longtin  et.  al.  1992];  such  a  distribution  has  also  been 
observed  experimentally. 
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Fig  8.  ISIHs  computed  via  numencal 
simulation  of  (9)  with 
(P ,  q  ,  0) ,  a,2 )  =  ( l  .6056 . 0.304 ,  ni  10 . 0. 1 34 1 
and  uncertainty  (see  text)  A;  =0 
(solid  curve),  0.1  (dotted  curve), 
and  0.25  (data  points).  Note  the 
transition  from  peaks  at  *772  to 
peaks  at  nT  for  increasing  At . 


Our  model  is  seen  to  reproduce  all  substantive  features  of  the  expcnmental  data:  in  addition 
to  the  characteristic  T-dependent  locations  of  the  successive  peaks,  the  modal  decay  rates  (except 
for  the  first  few  peaks)  are  exponential.  Analog  simulations  show  [Longtin,  Buisara  and  Mox> 
1991;  Longtin  et.  al.  1992]  that,  on  a  semi-logarithmic  scale,  the  decay  constant  is  proportional 
to  the  modulation  amplitude  6  for  fixed  noise  intensity  a2,  with  a  qualitatively  similar  relation 
ship  obtained  between  the  decay  rate  and  the  noise  strength  for  constant  8.  This  is  not  too  suppos¬ 
ing  since  the  noise  and  signal  arc  on  equal  footing  in  (9).  Wc  may  then  speculate  that,  over  a  cer 
tain  (as  yet  not  fully  defined)  range  of  parameters,  the  noise  and  signal  pla\  interchangeable  roie- 
in  determining  the  shape  of  the  ISIH.  Their  roles  are  not  completely  reciprocal,  howeser.  sik; 
the  peak-widths  in  the  ISIH  arc  dependent  on  a.\  Increasing  the  stimulus  amplitude  lead-  to  m 
increa.se  in  the  heights  of  the  lower  lying  peaks  This  is  consistent  with  experimental  observation 
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transcend  our  simple  description.  Clearly .  tiu>  is  a: 


which  meats  corisidcrat'ie  kiniiet  •  ’.;•!• 


Comparison  to  Integrate-l  ire  Models 

A  stochastic  model,  similar  in  spin!  to  the  deterministic  integrate-fire  model  sec  eg 
Keener.  Hoppcnstcadt  and  Rinzcl  19X1,  and  references  therein)  was  originally  developed  <  lei 
stein  and  Mandelbrot  1964|  to  try  to  explain  the  experimentally  observed  1SIH  corresponding  to 
spontaneous  firing  events;  as  pointed  out  above,  this  distribution  lunction  is  a  Gamma  distribu¬ 
tion  Assuming  the  underlying  dynamics  to  be  time-stationary,  a  random  walk  description  was 
invoked,  based  on  the  cornerstone  requirement  of  a  stable  distribution  function  for  the  probability 
density  of  first  passage  times  corresponding  to  the  dynamics.  The  state  variable  u-  was  assumed 
to  execute  a  biased  random  walk  to  an  absorbing  threshold  at  which  point  a  tiring  event  was 
designated  to  have  occurred  and  the  membrane  potential  u:  was  then  instantaneously  reset  to  its 
starting  value  (the  reset  mechanism  being  purely  deterministic  unlike  our  bistable  model,  m 
which  it  is  stochastic).  The  distance  between  die  origin  and  die  threshold  is  the  ''barrier  height"  z 
(analogous  to  the  height  U „  of  the  potendal  barrier  in  our  bistable  model)  in  die  Gerstein- 
Mandelbrot  description.  Further,  it  was  assumed  that  the  motion  in  phase  space  occurs  under  the 
influence  of  a  positive  drift  coefficient  p  which  was  defined  by  Gerstcin-Mandelbrot  as  the  differ¬ 
ence  between  die  drift  velocities  corresponding  to  excitatory  and  inhibitory  synaptic  inputs  (it  is 
neurophysiologicaily  reasonable  to  assume  these  velocities  to  be  different).  Then,  assuming  the 
presence  of  some  (as  yet  unquantified)  random  background  noise  which  is  taken  to  be  Gaussian 
delta-correlated  widi  zero  mean  and  variance  a2.  one  may  write  down  the  Langevin  equation  for 
this  process; 

u,  =  p  +  F(i),  (18) 


to  which  corresponds  the  Fokker  Planck  equation 


d  .  dP 
— P(U|.t).-n  — 


+  cr 


c )2P 
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(19) 


for  the  probability  density  function  P(u}  ,t).  This  equation  can  be  readily  solved  subject  to  the 
appropriate  boundary  conditions  and  the  probability  density  function  of  first  passage  times  writ¬ 
ten  down  in  the  form  [Gcrstein  and  Mandelbrot  1964], 
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(20) 


The  density  function  g(t )  reproduces  many  of  the  properties  of  experimentally  observed  ISIHs  for 
the  spontaneous  firing  case.  The  mean  first  passage  time  to  the  absorbing  threshold  is  calculated 
as  the  first  moment  of  g(t),  and  its  reciprocal  yields  an  average  firing  rate.  Variations  of  this 
model  incorporating  moving  boundaries  (which  mimic  refractoriness  and  are  therefore  closer  to 
neurophysiological  reality)  as  well  as  a  drift  term  that  is  linear  in  the  dependent  variable  u,  (the 
underlying  dynamics  is.  in  this  case,  representative  of  an  Omstein-Uhlenbeck  process),  have  been 
studied  by  Johanessma  [1968]  and.  Clay  and  Goel  [3973j. 


In  order  to  make  even  better  contact  with  experimental  results,  it  is  necessary'  to  provide 
reasonably  good  numerical  values  for  the  drift  coefficient  p,  the  "barrier  height"  z  and  the  back¬ 
ground  noise  variance  a2  in  the  above  model.  A  first  attempt  to  do  so  (while  simultaneously  pro¬ 
viding  a  test  of  the  goodness  of  fit  of  the  model  to  neurophysiological  data)  was  carried  out  by- 
Berger  ct  al  i  199()|.  They  earned  out  an  experiment  aimed  at  recording  the  inter-spike-interval 
distribution  from  extra-cellular  recordings  on  the  cat  visual  cortex.  Having  obtained  the  experi¬ 
mental  ISIHs.  they  were  able  to  compute  the  equivalent  model  quantities  p  and  z  via  the  mean 
and  standard  deviation  of  the  experimentally  obtained  ISIHs,  assuming  a  fixed  background  noise 
variance  rr  While  we  do  not  give  any  further  details  of  the  experiment,  it  is  noteworthy  that, 
.'me  these  self -consistent"  values  of  p, :  and  rr  were  substituted  into  the  first  passage  time  pro¬ 
bability  denstix  function,  an  excellent  fit  of  the  model  i20)  to  the  experimental  ISIHs  resulted  In 
.■■■rqccn;  publication,  Berger  and  Pribram  [1992;  extended  their  work  to  incorporate  the 


ihe  lulls  coupled  N-bod>  swem  ssithm  the  constraints  ol  the  theory  The  approach  lead"  n  - 
macroscopic  "potential  function  (•’  (defined  in  (loo  which  guarantees  gh  hal  stability  ol  uu: 
dynamic  system  lor  positive  a  without  having  to  constrain  ourselves  to  a  s.mmcirn  «.oupi;r.p 
matrix  J.  Our  theory  is  seen  to  agree  well  with  large  numerical  simulations  of  the  coupled  mo 
chastic  differential  equations  < 1  i.  within  the  hounds  imposed  by  the  constraints  < 2  ?  and  i  1 2  / 

The  theory  described  here  enables  us  to  describe  a  network  of  nonlinear  oscillator"  with 
nonlinear  coupling.  However,  because  of  the  separation  ol  time-scales,  the  dendritic  hath  t- 
tacitly  assumed  to  be  very  close  to  its  steady  state.  This  leads  to  the  quasi  hneari/anon  approxi  • 
relations  (9)  and  ill)  and  bnngs  our  description  of  the  bath  closer  to  other  quasi-lmear  dendnitc 
models  [see  e  g.  Segev  ct  al.  1989).  In  effect,  we  have  assumed  that  the  dendritic  patches  are 
only  weakly  bistable:  they  are  not  "strong"  threshold  devices.  These  assumptions  also  bring  our 
approach  closer  to  conventional  mean-held  theories  [see  c.g.  Amit  1989  for  an  overview)  1‘nlikc 
such  theories,  the  current  approach  does  not  depend  on  a  large  number  A  of  entities  to  improve 
its  convergence  although,  as  pointed  out  earlier,  in  the  presence  of  additional  elements  in  < 1  >  w  ith 
similar  time  constants,  recourse  to  a  more  conventional  mean-field  approach  may  be  unavoidable 
Further,  it  is  interesting  to  note  that  the  notion  of  representing  the  dendritic  bath  as  a  tcsselation 
of  elemental  volumes  each  desenbed  by  a  quasilinear  stochastic  differential  equation  for  an 
activation  function  u ,  (i>l)  is  similar  in  spirit  to  existing  compartmcntal  models  of  dendntte  trees 
(see  e.g.  Segev  et.  al.  1989).  Our  results  appear  to  be  independent  of  the  choice  of  the  statistics  of 
the  elements  of  J;  repeating  the  calculations  of  this  paper  with  the  J,,  drawn  from  a  uniform  dis¬ 
tribution  yields  qualitatively  similar  results  although,  as  pointed  out  earlier,  Gaussian  statistics 
may  be  more  reasonable  from  a  neurophysiological  perspective.  In  this  connection  it  is  worth 
pointing  out  that  one  expects  typically  small/sparse  interactions  (charactenzcd  by  coefficients 
Jq  —>OiJ>  1)  betw'een  dendritic  volumes  so  that  these  coefficients  may  indeed  be  reasonably 
characterized  by  a  sharply  peaked  (about  zero  mean)  Gaussian.  The  distribution  of  the 
coefficients  JX/  and  (these  coefficients  characterize  the  interaction  between  the  soma  and  the 
dendritic  volumes)  is  broader.  Good  agreement  between  the  probability  density  P(ux)  for  the 
reduced  system  (9)  and  the  exact  system  (1)  is  also  obtained  for  somewhat  larger  q  values  and 
noise  strengths  (within  the  bounds  of  the  inequality  n2)),  although  the  agreement  begins  to  break 
down  when  the  adiabatic  condition  on  the  frequency  is  violated.  The  magnitudes  and  signs  of  the 
J,,  can  be  very  important  in  determining  the  overall  sign  of  the  renormalized  coefficient  [3  in 
(lib)  and  this,  in  turn,  determines  the  modality  of  the  potential  (10).  For  a  monosiable  potential 
(P  <  a)  the  cell  body  is  always  quiescent  and  there  are  no  cooperative  effects.  The  bath  coupling 
can  render  a  monoscabfe  potential  (for  the  isolated  cell  body)  bistable,  under  certain  conditions, 
thereby  imparting  a  firing  capability  to  the  neuron;  the  opposite  effect  can  also  occur.  This  is  evi¬ 
dent  from  the  definitions  (llb.c):  changing  the  dendritic  parameters  changes  the  barrier  height 
and  the  location  of  the  elliptic  points  of  the  effective  potential  (10)  that  characterizes  the  soma 
dynamics,  while  also  renormalizing  the  modulation  amplitude.  These  changes  lead,  in  turn,  to 
changes  in  the  SNR  given  by  the  expressions  (16)  and  (17). 

The  approach  to  the  processing  of  information  in  noisy  nonlinear  dynamical  systems,  based 
on  the  probability  density  of  residence  times  in  one  of  the  stable  states  of  the  potential  offers  an 
alternative  to  the  FFT,  and  has  been  applied  [Longtin,  Bulsara  and  Moss  1991;  Longtin  ct  al 
1992)  in  the  theoretical  construction  of  inter-spike -interval  histograms  (ISIHs)  that  describe  neu¬ 
ronal  spike  trains  in  the  central  nervous  system.  This  model  exhibits  remarkable  agreement  with 
data  obtained  in  two  different  experiments  some  25  years  apart  [Rose  et.  al.  1967;  Stegal  1990)  as 
well  as  with  the  more  recent  data  of  Rhode  [1991;  unpublished];  figures  6  and  7  demonstrate  this 
agreement.  The  approach  of  Longtin  et  al.  has  been  contrasted  with  more  conventional  theories 
of  ISIHs  based  on  integrate-and-ftrc  (IF)  models  in  which  the  activation  performs  a  random  walk 
to  an  absorbing  barrier  and  is  then  reset  to  its  initial  value.  In  the  absence  of  an  absolute  refrac¬ 
tory  penod,  the  two  approaches  may.  in  fact,  converge  with  the  mean  firing  rate  (computed  j"  the 
reciprocal  of  the  mean  first  passage  time)  in  the  IF  model  corresponding,  roughly,  to  the  mean 
duration  of  a  full-cycle  switching  event  in  the  bistable  diffusion  model  of  l  ongtin  et  a!  '  he 
approach  of  Longtin  et  al  .  however,  seems  to  offer  the  most  elegant  treatment  of  the  [SfH.  cer 
minis  it  permits  one  to  match  the  model  with  experimental  data  (far  more  closels  than 
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